            real fpair = dwdr + dudr;

            vecscaleaddv(icell->f[i], icell->f[i], datom, 1, fpair);
            vecscaleaddv(jcell->f[j], jcell->f[j], datom, 1, -fpair);
            vec<real> f;
            vecscale(f, datom, fpair);

            stat->evdwl += w;
            stat->ecoul += u;
          }
        }
      }
    }
  }